Morales-Castilla et al. (2015, TREE) Inferring biotic interactions from proxies.
Vanessa, Gian-Marco, Thomas
June 16, 2015
Introduction
This is the reproduction of the analyses presented in the article Inferring biotic interactions from proxies by Ignacio Morales-Castilla, Miguel G. Matias, Dominique Gravel, and Miguel B. Ara??jo (2015, TREE). The article can be found here.
Load data
Clear R environment and load libraries.
# Clear workspace.
rm(list=ls())
# Load libraries.
library(igraph)
library(RCurl)## Loading required package: bitopsLoad the data from GitHub straight into R.
Serengeti food web.
# Species in the Serengeti food web.
dd <- getURL("https://raw.githubusercontent.com/opetchey/RREEBES/master/MORALES-CASTILLA_etal_2015_TREE/Data/Serengeti%20data/journal.pcbi.1002321.s004.CSV",
cainfo = system.file("CurlSSL", "cacert.pem", package = "RCurl"))
dd.01 <- read.csv(text = dd)
# Feeding links in the Serengeti food web.
dd <- getURL("https://raw.githubusercontent.com/opetchey/RREEBES/master/MORALES-CASTILLA_etal_2015_TREE/Data/Serengeti%20data/journal.pcbi.1002321.s005.CSV",
cainfo = system.file("CurlSSL", "cacert.pem", package = "RCurl"))
dd.02 <- read.csv(text = dd)
# Consensus partition.
dd <- getURL("https://raw.githubusercontent.com/opetchey/RREEBES/master/MORALES-CASTILLA_etal_2015_TREE/Data/Serengeti%20data/journal.pcbi.1002321.s006.CSV",
cainfo = system.file("CurlSSL", "cacert.pem", package = "RCurl"))
dd.03 <- read.csv(text = dd)
# Link densities between groups in the consensus partition.
dd <- getURL("https://raw.githubusercontent.com/opetchey/RREEBES/master/MORALES-CASTILLA_etal_2015_TREE/Data/Serengeti%20data/journal.pcbi.1002321.s007.CSV",
cainfo = system.file("CurlSSL", "cacert.pem", package = "RCurl"))
dd.04 <- read.csv(text = dd)Calculate the number of species and the number of potential links one could expect from this.
# Species only in Cuban food web.
n.species.Serengeti <- length(dd.01[,1])
n.links.species.Serengeti <- n.species.Serengeti^2-n.species.SerengetiNote:
- There are 175 species in the Serengeti food web. This results in 25760 potential links (\(n^2-n\)).
The cuban coral reef food web.
dd <- getURL("https://raw.githubusercontent.com/opetchey/RREEBES/Morales-Castilla_etal_2015_TREE/MORALES-CASTILLA_etal_2015_TREE/Data/Coral%20reef%20data/857470.item.4%20-%20Guild%20key%2C%20all%20data%20sets.csv",
cainfo = system.file("CurlSSL", "cacert.pem", package = "RCurl"))
dd.05 <- read.csv(text = dd)
str(dd.05)## 'data.frame': 755 obs. of 6 variables:
## $ Guild.Number : int 1 2 2 2 2 2 2 2 2 2 ...
## $ Taxa : Factor w/ 753 levels ""," ",..: 1 661 120 126 127 132 282 396 407 685 ...
## $ Fish.Body.Length: num NA NA NA NA NA NA NA NA NA NA ...
## $ Cayman.Islands : Factor w/ 3 levels ""," ","x ": 1 1 1 1 1 1 1 1 1 1 ...
## $ Cuba : Factor w/ 3 levels ""," ","x": 1 1 1 1 1 1 1 1 1 1 ...
## $ Jamaica : Factor w/ 3 levels ""," ","x": 1 1 1 1 1 1 1 1 1 1 ...# Convert 'dd.05$Guild.Number' frmo integer to factor.
dd.05$Guild.Number <- as.factor(dd.05$Guild.Number)
# str(dd.05)
# Only data with species from the Cuban food web.
# Guilds only in Cuban food web.
dd.05.a <- subset(dd.05, Cuba=="x")
n.guilds.Cuba <- length(levels(factor(dd.05.a$Guild.Number)))
# Species only in Cuban food web.
n.species.Cuba <- length(subset(dd.05, Cuba=="x")$Cuba)
n.links.species.Cuba <- n.species.Cuba^2-n.species.Cuba
# Species and guilds in all three food webs.
dd.05.b <- subset(dd.05, Cuba=="x" | Jamaica=="x" | Cayman.Islands=="x")
n.species.all <- length(dd.05.b[,1])
n.links.species.all <- n.species.all^2-n.species.all
n.guilds.all <- length(levels(factor(dd.05.b$Guild.Number)))
n.links.guilds.all <- n.guilds.all^2-n.guilds.allNote:
Well, overall there are 175 species in the three food webs, and 171 species in the Cuban food web. This would give a total number of potential links of 30.450 or 29.070, respectively.
There are 144 guilds in all three food webs and 140 guilds in the Cuban one. Based on guilds, this would give a a total number of links of 20.592 or 19.460, respectively. The question is: where do the 231 species/guilds come from that Ingacio Morales-Castilla was telling us entered their calculations…?!?
How to…?
“The procedure for building the backbone of an interaction network starts with the identification of species more likely to share similar interactions. The concept is similar to that of modules [28], but we avoid this terminology because modules are usually determined a posteriori and can also refer to simple assemblages of species such as linear food chains or apparent competition [29]. Instead, we define interacting groups based on a priori expectations of interactions. The concept is also analogous to that of guilds [30]. Guilds, however, are restricted to species shar- ing similar resources, and thus do not encompass non- consumptive interactions such as competition or niche construction [27]. A flexible definition of interacting groups based on traits, phylogenies, and geographical distribu- tions would enable combination of heterogeneous information.” (page 4, bottom-right)
“Interacting groups of species are defined a priori to simplify the removal of forbidden links. The groups were defined based on the trophic hierarchy of the different species within each eco- system (e.g., primary producers, grazers, small and large carnivores). This process of trophic classification of species led to identification of forbidden links and removal of 30% of all potential direct links in the coral reef, and 22% in the Serengeti (e.g., herbivores eating predators; Figure 1).” (page 5, bottom-left)
“Refinement of the species groupings was achieved by con- sidering the characteristics of the consumer species (e.g., distinguishing small versus large carnivores in the Seren- geti example, or separating invertebrate feeders, omnivo- rous, and carnivorous fish in the coral reef example; Figure 1).” (page 5, bottom-right)
Serengeti foob web analysis
Taking all potential links into account the food webs would look like this:
dd.01.a <- expand.grid(dd.01$code, dd.01$code)
# Build the graph.
f <- graph.data.frame(dd.01.a, directed=F, vertices=NULL)
# Now, remove inraspecific trophic interactions - there is no cannibalism.
dd.01.b <- dd.01.a[which(dd.01.a[,1]!=dd.01.a[,2]),]# Now, plot the graph.
plot.igraph(f,
layout=layout.circle,
main="Serengeti food web (no information)",
vertex.size=8,
vertex.color=adjustcolor("#A6B07E", alpha.f=0.6),
vertex.label.cex=0.5,
vertex.label.family = "sans",
edge.width=0.2,
edge.color=adjustcolor("#999999", alpha.f=0.3),
vertex.label.color="#7D1D3F",
vertex.frame.color="white")
This is what the real food web looks like. So far, there is no ordering according to any traits (e.g., nutrition, size…)
# Build the graph.
g <- graph.data.frame(dd.02, directed=F, vertices=NULL)
# Display the food web graph.
# tkplot(g)
# plot.igraph(g)plot.igraph(g,
# layout=layout.fruchterman.reingold,
layout=layout.circle,
main="Serengeti food web (real food web)",
vertex.size=8,
vertex.color=adjustcolor("#A6B07E", alpha.f=0.6),
vertex.label.cex=0.5,
vertex.label.family = "sans",
edge.width=0.2,
edge.color=adjustcolor("#333333", alpha.f=0.6),
vertex.label.color="#7D1D3F",
vertex.frame.color="white")
Once we know about the trophic levels the species/guilds belong to we can use layout.reingold.tilford(g, flip.y=FALSE) to introduce some hierarchy. For example, if we assign three trophic levels just randomly then the food web can be displayed like this:
plot.igraph(g,
# layout=layout.fruchterman.reingold,
layout=layout.reingold.tilford(g, flip.y=FALSE),
main="Serengeti food web (no trait information)",
vertex.size=8,
vertex.color=adjustcolor("#A6B07E", alpha.f=0.6),
vertex.label.cex=0.5,
vertex.label.family = "sans",
edge.width=2,
vertex.label.color="#7D1D3F",
vertex.frame.color="white")